"Cooling by heating" - demonstrating the significance of the longitudinal specific heat 



(N 

o 

(N 



Jon J. Papini, Jeppe C. Dyre, and Tage ChristensenEI 

DNRF Centre "Glass and Time ", IMFUFA, Department of Sciences, 
Roskilde University, Postbox 260, DK-4000 Roskilde, Denmark 
(Dated: January 29, 2013) 

Heating a solid sphere at the surface induces mechanical stresses inside the sphere. If a finite amount of 
heat is supplied, the stresses gradually disappear as temperature becomes homogeneous throughout the sphere. 
We show that before this happens, there is a temporary lowering of pressure and density in the interior of the 
sphere, inducing a transient lowering of the temperature here. For ordinary solids this effect is small because 
c p — cv ■ For fluent liquids the effect is negligible because their dynamic shear modulus vanishes. For a liquid 
at its glass transition, however, the effect is generally considerably larger than in solids. This paper presents 
analytical solutions of the relevant coupled fhermoviscoelastic equations. In general, there is a difference be- 
tween the isobaric specific heat, c p , measured at constant isotropic pressure and the longitudinal specific heat, 
C(, pertaining to mechanical boundary conditions that confine the associated expansion to be longitudinal. In 
the exact treatment of heat propagation the heat diffusion constant contains ci rather than c p . We show that the 
key parameter controlling the magnitude of the "cooling-by-heating" effect is the relative difference between 
these two specific heats. For a typical glass-forming liquid, when temperature at the surface is increased by 1 K, 
• a lowering of the temperature in the sphere center of order 5 mKis expected if the experiment is performed at 

CN| ' the glass transition. The cooling-by-heating effect is confirmed by measurements on a 19 mm diameter glucose 

sphere at the glass transition. 

^ PACS numbers: 

O 

I. INTRODUCTION 

Most solids and liquids expand when heated. Heat diffusion is a notoriously slow process, and heating a solid sample at its 
surface induces stresses in the sample that only disappear when temperature gradually becomes again homogeneous throughout. 
Heating a lightly fluent fluid that has a free surface (i.e., is free to expand), on the other hand, makes the entire sample expand on 
the time scale set by the sound velocity and sample dimensions. In this case there are no transient stresses beyond the acoustic 
time scale. A liquid close to its glass transition provides an interesting case in between solid and fluid behavior. Such a liquid 
■ behaves like a solid on time scales shorter than the Maxwell relaxation time tm = f//Goo where r/ is the shear viscosity and Goo 
the instantaneous shear modulus. The Maxwell relaxation time becomes longer than one second when a liquid approaches its 
calorimetric glass transition, implying that induced stresses survive for seconds or more. 

This paper discusses the "cooling-by-heating" effect that arises when a sample is heated at a free surface. We show that this 
effect, which is present in all hard solids with a non-zero thermal expansion coefficient, is generally magnified considerably for 
glass-forming liquids close to their glass-transition temperature T g . This is because close to T g the liquid is solid-like by having 
a large, non-zero dynamic shear modulus on short time scales and, at the same time, is liquid-like by having a fairly large thermal 
expansion coefficient. 

Returning to the case of a solid, what happens when heat is supplied at the (free) surface of a spherical sample? The outermost 
layers attempt to expand, obviously, but a priori one may imagine two different possibilities: 1) the expansion presses inwards, 
resulting in an increase of the pressure at the center of the sphere; or: 2) the expansion turns outwards, thus transmitting a 
negative pressure into the sphere. Which of the two possibilities that applies is answered by the application of standard thermo- 
elasticity theory to the problem of calculating the stresses induced by the heating. This is done in the present paper. It turns out 
that case 2 applies — the sphere expands and pressure decreases in the interior of the sphere. This induces an adiabatic cooling 
inside the sphere. The phenomenon of cooling caused by heating at the surface is referred to below as the cooling-by-heating 
effect. 

The solution of the coupled thermomechanical equations detailed in Sec. [HI] shows that the cooling-by-heating effect is 
proportional to the difference between the reciprocals of the isobaric specific heat, c p and the longitudinal specific heat, q (all 
specific heats are per unit volume); the latter quantity was introduced and discussed in Refs. QHH The longitudinal specific heat 
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is related to the isochoric specific heat, cy, by 



where A/5 and Mt are the adiabatic and isothermal longitudinal moduli respectively. This is analogous to the standard thermo- 
dynamic relation c p = [Kg / Kt)cv relating the isobaric specific heat to the isochoric specific heats in terms of the adiabatic, 
Ks and isothermal, Kt bulk moduli. Since Ms = K$ + (4/3)G and Mt = Kt + (4/3)6?, where G is shear modulus, one 
readily finds that c; is in between cy and c p . As we shall see, the relative difference a; = (c p — c\)j c p controls the strength of the 
cooling-by-heating effect, and we thus term this quantity the "longitudinal thermomechanical coupling constant". Combining 
the equations above ai is found to be the product of two factors JU, 

Cp — ci 4 G Cp - cy ... 

ai = — = - — — ■ (2) 

c p 6 Mt c p 

The first factor shows that there is only cooling by heating if the shear modulus is non-vanishing compared to the longitudinal 
modulus. The other factor, "the thermomechanical coupling", a = (c p — cv)/c p is a dimensionless measure of the coupling 
between thermal and mechanical perturbations. It can be expressed in terms of the expansivity, a p = (l/V)(dV/dT) p , as 
follows: 

c P - c v Tqo?K t 
a = -z = ^ , (3) 

Cp Cp 

where To is the temperature. It follows that the cooling -by-heating effect is quadratic in the thermal expansion coefficient a p . 

Since solids typically expand significantly less upon heating than do liquids, the cooling-by-heating effect is generally small 
in solids. As an example, for solid glucose the thermal expansion coefficient |3] is 1.1 • 1CP 4 K _1 close to the glass transition 
whereas for liquid glucose it is 3.7 • 10~ 4 K _1 in the same temperature region. This potentially enhances the cooling-by-heating 
effect by a factor of 11. However the changes in c p from 1.91 • 10 6 JK _1 m~ 3 to 3.05 • 10 6 JK _1 m -3 and in Kt 0] from 
10.75 • 10 9 Pa to 6.49 • 10 9 Pa reduces this to a factor of 8. Here we have used a density of 1.52 • 10 3 kg m~ 3 to convert specific 
heat data from mass to volume. It is, however, not unusual for liquid expansivities to be near 10 -3 K _1 for which we would 
expect an enhancement of the thermomechanical coupling a = (c p ~ cv)/c p by a factor of 30. The shear modulus of glucose in 
the glassy state is G^ = 3.1 • 10 9 Pa as deduced from the shear compliance data of Meyer and Ferry 

The above relations all generalize to deal with complex, frequency-dependent (dynamic) specific heats and moduli and expan- 
sivity, which are the relevant quantities when studying glass-forming liquids. Near the glass transition the cooling-by-heating 
effect may be studied on second time scales. Here, upon increasing the frequency, the factor G/Mt of Eq. (O increases while at 
the same time the factor (c p — c v )/c p decreases. The enhancement of the cooling-by-heating effect is thus critically dependent 
on the relative time scales of the different relaxation processes at the glass transition. If the shear stress relaxes faster than the 
the volume processes, the cooling-by-heating effect may not be pronounced. This situation is illustrated in Fig. Q] The model 
describing the relaxation behavior between high- and low-frequency values is described in section [IV] 

The present work discusses the basis of cooling by heating by referring to the equations of standard linear thermoviscoelastic- 
ity. Section|II]introduces the general framework of thermoelasticity and thermoviscoelasticity. It is shown that the heat diffusion 
constant involves the longitudinal specific heat. Section[III]discusses the case when a finite amount of heat is fed into a sample 
at its surface at t = 0, as well as the experimentally easier realized case when temperature is suddenly increased at the surface. 
That section also presents analytical calculations of the ordinary solid case for which the constitutive properties do not undergo 
relaxation. Section [IV] gives calculations of a model glass-forming liquid, i.e., the case when the constitutive properties are 
frequency dependent. We estimate the effect to be of order 5 mK in the center of a sphere for a temperature increase at the 
sphere surface of one Kelvin. Section [Vl confirms this prediction for measurements on a glucose sphere. Sections VI and IVIII 
briefly discuss and summarize the paper. 



II. THERMOELASTICITY AND HEAT DIFFUSION. 



Thermoelasticity deals with problems where displacement field u(r, t) and temperature field T(r, t) couple. It is a linear theory 
of small deformations given in terms of the strain tensor e ^ = | ( ^f- + ) and small temperature increments ST = T — Tq 
relative to a reference temperature Tq. The material properties of a thermoelastic medium is given by the linear constitutive 
equations that expresses stress and increments in entropy density <5^in terms of and ST. The hydrostatic pressure is related 
to the trace of the stress tensor p = — 1/3 an and the relative compression is the trace of the strain tensor e = J^i e u = V • u. 
The following constitutive equations lg] define the shear modulus G, the isothermal bulk modulus, Kt, the isochoric specific 
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FIG. 1: Sketch of overlapping relaxations of AG/ (3Mt) (green) and the thermomechanical coupling, a — Ta p Kr /c p (red). The longitudinal 
coupling constant a; = (c; — c p )/c p (black) is the product of those two frequency dependent functions. Near the crossing of the two curves 
the difference between c; and c p is generally largest. This determines the time scale of experiments on glass-forming liquids, where the 
cooling-by-heating effect is particularly large. 



heat, cy and the isochoric pressure coefficient /3y: 



crij +pSij = 2G{e l3 - ^eSij) (4) 
p = -K T e + f3 v 5T (5) 
5j=(3 v e+^6T. (6) 

We follow Biot J3l in assigning the symbol j3 to the thermodynamic pressure coefficient 

The material is furthermore characterized by the heat conductivity, A, which enters Fourier's law for the entropy current density 

is- 

h=-^VT. (8) 

The interest in thermoelastic problems has since Duhamel I8J mostly been focused on the calculation of thermal stresses 
deriving from an evolving temperature field. In the classical thermoelasticity theory the displacement field and temperature fields 
are partially decoupled [|9lll0|]. This comes from assuming that the development of the temperature can be found independently 
of the stresses by the conventional heat-diffusion equation: 

dST 2r 

— = DV 2 5T . (9) 
at 

Here D is a heat diffusion constant. After solving this equation the displacement field can be found from the quasi-static stress 
equilibrium equation: 

M T V(V-u) -GV x (V xu) -f3 v V5T = (10) 
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This approximate theory is referred to as the theory of thermal stresses OJ. According to many authors the correct 

treatment appeared remarkably late in the development of thermoelastic theory with Biot's paper 101 in 1956. Lessen lll4"tl 
considered similar problems the same year. The heat diffusion equation Eq. © is replaced by 

dST dV-u 2 

cv-xr+T Pv — k— = AV<5T, (11) 
at at 

which follows from entropy conservation 

s~ v '* <12) 



when this is combined with Eqs. (|6]l and ([8). Entropy conservation may seem strange at first sight, but the entropy production 
per volume associated with heat conduction is — j • -r- VT = 4(VT) 2 , i.e., a second-order effect disappearing in a linearized 

1 1 

theory. 

In most cases the ordinary, decoupled heat-diffusion equation is a good approximation in the manner it is used in the theory 
of thermal stresses. However, this approximate theory is not able to describe the phenomenon of cooling by heating, which is 
the theme of this paper. It should be noted that the heat-diffusion equation with the diffusion constant containing the isobaric 
specific heat c p is exact for the non-viscous liquid state or soft matter with G = 0, if part of the boundary (with normal vector 
n) is free to expand, i.e., . a^rtj = 0. The proof runs as follows: The assumption G/Mt = simplifies Eq. ( TTOb to 

V(K T V -u- f3 v ST) =0. (13) 

However, the terms under the gradient is according to Eq. (|5} nothing but minus the pressure increment. Thus we conclude this 
pressure increment is uniform in space and only depends on time. Moreover, Eq. (0) ensures that all diagonal elements of the 
stress tensor are identical and equal to minus this pressure increment. If the normal component J2j a ij n j is zero on P art °f me 
boundary, it follows that the pressure is also zero there, but then it is zero throughout the body. Equation ((5) then reduces to 
V • u = a p 5T. Inserting this in the entropy equation Eq. (fTTT i. one arrives at the ordinary decoupled heat diffusion equation 
with D, p = X/c p 

O frri 

— =D p V 2 5T, (14) 

when noticing that c p — c v + T j^ = cy + T a ap\K T . 

As se have seen, the temperature field in general does not exactly obey a diffusion equation. It does so when c p = cy 
fly = 0) or for certain boundary conditions when G = 0. However, as emphasized by Biot J3l, the entropy density in 
fact does fulfill a diffusion equation and moreover with a diffusion constant containing the ubiquitous longitudinal specific heat: 
Applying the divergence operator to the inertia-free stress equilibrium equation Eq. (TTOb . gives 

Vh = ^fv 2 8T. (15) 

Applying the Laplacian to the constitutive equation © yields 

V 2 ^=(|k + ^)v 2 5T. (16) 



Fouriers law and the entropy conservation Eq. (TTZt gives 



9 4=-VHT, (17) 
at T 



and thus 

£=-V a * (18) 
at ci 

with q = cy + Tojj^ being the longitudinal specific heat 0101. 

Note that this result is limited to the inertia-free cases. If one wishes to study coupled mechanical and thermal waves, the 
inertia-term P^-u must be added on the right side of Eq. ( fTOt . Solutions of the equations in this case have been studied 
extensively (gy also in the spherically symmetric situation. Note, however, that acoustic wavelengths are much longer than 
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thermal wavelengths. Thus for a sample of a certain size there is an interesting time regime where acoustic waves have settled, 
but thermal diffusion has barely begun. Take as an example a sphere of radius 1cm. For a typical sound velocity of 10 3 m/s and 
heat diffusion constant of 10~ 7 m 2 /s, the sound traveling time is 10 fis while the diffusion time is 1 ks. It is within this time 
regime we will find the cooling-by-heating phenomenon. Although the solution restricted to the inertia-free case that we present 
below is in principle contained in the coupled acoustic-thermal wave solutions including inertia, the phenomenon is obscured by 
the complicated structure of these solutions and seems not to have been recognized. 

The thermoelastic theory that was originally developed for elastic solids without relaxation is easily extended to a thermovis- 
coelastic theory taking relaxation of all the constitutive parameters into account, as it is necessary for relaxing liquids near the 
glass transition. The most straightforward way of generalizing is to interpret the equations in the frequency domain allowing all 
the constitutive parameters to be complex functions of the angular frequency to. The cases we study in the frequency domain 
cover thus both solids and thermoviscoelastic liquids, but can only be transformed analytically into the time domain for solids. 
For relaxing liquids one must do the transformation numerically. 



III. ANALYTICAL SOLUTIONS OF THE SPHERE-HEATING PROBLEM 



A. The case when the heat flow is controlled at a mechanically free boundary 

This subsection presents the analytical solution in the frequency domain to the situation when a sphere of a general viscoelastic 
material is subjected to a periodic heat input at the surface. The solution shows the temperature in the center at high frequencies 
varying 180° out-of-phase with respect to the heat oscillation at the surface, indicating the cooling-by-heating effect. In order to 
give a more lucid and transparent understanding of the phenomenon we translate the solution to the time domain. This can be 
done analytically by an inverse Laplace transformation if there is no frequency dependence of the constitutive properties. That 
is, we calculate the temperature and stress profile throughout the sphere following a heat-step input at the surface at time zero. 

Consider the case when a periodically varying heat SQ(t) = Re {5Qe lult } is supplied at the surface of a sphere of ra- 
dius R. The surface is assumed to be mechanically non-clamped, i.e., the sphere is free to expand. This translates into the 
boundary condition that the radial component of the stress tensor is zero at the surface, a rr (R,oj) = 0. We wish to calculate 
how the periodically varying temperature and displacement fields vary throughout the sphere, i.e., to calculate the complex 
frequency-dependent amplitudes of temperature, 8T(r,u>), and radial displacement field, u(r, oj). From these quantities the 
stress components oy r (r, oj), etc, may be calculated. 

Denoting the angular frequency by oj, the position vector by r, the complex frequency-dependent radial displacement field by 
u(r, oj) = u(r, oj)r/r, the coupled thermoelastic equations ( TTOb and (fTTl i become 



|- M T r- 2 -^-(r 2 u)-(3 v ST 
d 

(iuj)cvST + (ioj)T /3v r~ 2 — (r 2 u) 

or 



= 



(19) 
(20) 



The four boundary conditions are: 

1. No displacement at the center: u(0, oj) = 0; 

2. No temperature gradient at the center: ^p-(0, oj) = 0; 

3. Free surface, i.e., no radial stresses at the surface: a rr (R, oj) = 0; 

4. Heat supply boundary condition at the surface: A^r(R, oj) = ioj ■ 
Denote the volume of the sphere by Vq 



yjiojci (oj)/X. Define furthermore the functions 

h{r/R,k 2 R 2 ) 
f 2 (r/R,k 2 R 2 ) 

Introduce the characteristic heat diffusion time 



\ttR and define the complex frequency-dependent thermal wavevector by fc 

(fci?) 3 



1 sinh(fcr) 
3 kr kRcosh(kR) — smh(kR) 
R^ 3 (kr) cosh(fcr) — sinh(fcr) 



(fci?) cosh(fci?) - sinh(fci?) 



(21) 
(22) 



r = (c 4 (oj)/X)R 2 . 



(23) 
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Then one has k 2 R 2 



iuir and the solutions to Eqs. ( fT9l ) and d20b are 

1 



5T(r,u>) 
u(r, uj) 



V ci 
1 a 



{-ai + fi(r/R,iwT)}5Q 
4 G K s 



f 2 {r/R,iwT)\5Q 



(24) 



(25) 



3V c p [3M S M s 

These solutions were found by the transfer-matrix approach (see Ref. @ and Appendix IA j} , and can be verified by insertion 



noticing that fx (p, s) 



°(p*f 2 (p,s)), £fx(p,s) = hsf 2 (p, S ) and X §-o 2 =sfx(p,s). 



i 

3P 1 dp V -^P) "))' dpJ 1 ^' °> ~ zr°Jt\Hi "J -pi Qph 1 QpJ 

Consider the low- and high-frequency limits of these expressions. The functions fx and f 2 both have the limits 1 for u —> 
and for ui — > oo. Thus, as expected ST — > SQ/(VqC p ) in the low-frequency limit when heat has had time to distribute 
throughout the sphere. At high frequency the temperature amplitude becomes 



ST(r) 



1 



V ci 



ai5Q 



1(1 

V y Cl 



1 



)SQ for 



(26) 



If we for a moment consider the non-relaxing case where the specific heats are real, we see that the temperature amplitude is in 
counter-phase to the heat amplitude since q < c p . For a propagating thermal wave it would not be surprising that temperatures 
at some distance - e.g., at a half wavelength - had opposite phases. However, Eq. ( |26] i holds throughout the sphere and is not 
associated with the diffusive thermal wave. This will be even more clear when we consider the response to a heat step input later 
on. 

We see that the longitudinal coupling constant a; controls the magnitude of the cooling-by-heating effect. The ratio of the 
amplitudes of the temperature in the center and the heat input at the surface is shown in Fig. [2] The phenomenon "cooling by 
heating" is indicated at high frequencies, albeit this is more conspicuous in the time domain. 




1 2 

log 10 (a)T) 



FIG. 2: The real part of the ratio between the complex amplitudes of temperature at the center and the heat supplied at the surface scaled with 
the isobaric heat capacity. At high frequencies the limit becomes the negative value — ai/(l — ai) — 1 — c p /c; . 



For the displacement field we find for the low-frequency limit the natural result, 



. , 1 SQ 
u[r) -> ^a p —- 
3 V c p 



r for uj 



(27) 



determined by the final temperature rise and the linear thermal expansion \a p . However, at high frequencies we find 



, 1 5Q AG 



(28) 



Notice that this displacement, which is responsible for the cooling-by-heating effect, is only present when G ^ 0. 

In the spherically symmetric case there are only two different components of the stress tensor, a rr and <jqq = cr^. It follows 



from Eqs. © and © andjhe fact that e rr = and e ee = = f that a rr = {K T + + 2(K T - §G)^ - (3 V 8T, 



which by Eqs. d24l i and d25l l becomes 



a rr (r,uj) = ^ V {1 - f 2 (r/R,iujT)}5Q . 

3 Mt Vqci 



(29) 
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Likewise, a gg = (K T - § + 2(iT T + \G)^ - /3yc5T, which becomes 

v 6g {r,u,) = \-£ r p-{l- \h(r/R,iur) + \h{r/R, iur)\ SQ . (30) 
3 Mt Vqci I 2 2 J 

Thus at high frequencies there is an isotropic, uniform tensile stress in the interior of the sphere of the magnitude 

<7 rr , V88 , O>0 -> TF~^Q for 0J -> OO (31) 

3 iWr VoQ 

On the other hand, all stresses vanish for w — > (as expected). 

In order to gain insight into the cooling-by-heating effect and show the significance of the longitudinal coupling constant ai 
we transform Eqs. (l24l . ( |29l l and ( |30l into the time domain, but only for a solid, i.e., in the case when all constitutive properties 
are frequency independent. If a delta function heat flux is applied at t = 0, the heat supplied at the surface is a Heaviside step 
function, 5Q(R, t) = 8QoH(t); in this case calculating the inverse Laplace-Stieltjes transform leads to the following expressions 
for the temperature and stresses as functions of time after t = (see Appendix IA 21 ): 

5T(r,t) = ^-{- aj + F 1 (r/i?,t/r)}5Qo, (32) 

°Vr(r,t) = \-^-^{l-F 2 {r/R,t/T)}5Q^ (33) 
3 M T V ci 

*ee(r,t) = ~p- {l - ^(r/^i/r) + \f 2 {t/R, t/r)\ 6Q , (34) 
3 Mt VoQ 12 2 J 



where 



WiWr) = l + ^£!^e-^ (35) 
3 r f sm x n 

n— 1 v y 



F 2 (r/i?,Vr) = 1 + 2(^1 E ^(i^)-ix n cos(ix n ) e _ 4t/T _ 
V 7 '/ ^ ^sm(x„) 

Here x\ < x 2 < ■■■ are the positive roots of the transcendental equation x = tan(x). Note that F\ and F 2 — > for t — > 
and Fi and F 2 — > 1 for f — > oo. When a; = there is no cooling -by-heating effect according to Eq. (l32t . Furthermore, a; = 
implies either G = or (3y = and there is no immediate expansion and no induced stresses. 

When ai ^ the situation is quite different. In Fig. [3]we plot the scaled temperature change (c p Vo/5Q)ST(t/T;r / R) for 
a several radii r/R as given in Eq. d32b . The longitudinal coupling constant is here fixed to a; = 0.091 and time is given in 
units of the characteristic heat diffusion time r. The figure clearly shows the cooling-by-heating effect. Since a finite amount 
of heat was added at the surface at t = 0, the surface temperature initially diverges. The interior of the sphere, even close to 
the surface, instantaneously cools to a uniform temperature. The expansion of the surface is immediately felt in the interior, and 
since no heat has yet arrived by diffusion, it cools adiabatically. This initial response is followed by an evolution in time where 
the temperatures of the different parts of the sphere converge and eventually equilibrate. 

In order to understand better the physics of cooling-by-heating we consider the components of the stresses given by Eqs. (|33l 
and (|34i l, respectively. In Fig. |4]the a rr component of the stress tensor is plotted scaled with the initial uniform interior stress 
ctq = li^TyTj^Qo- First we note that the boundary condition a rr (R,t) — is fulfilled. As the surface receives heat and 
expands, an immediate traction is felt in the interior of the sphere. a rr is positive, seeking to stretch a volume element in the 
radial direction under the entire evolution to thermal equilibrium. 

The scaled stress component age (r, t), is shown in Fig. [5] One notices an immediate, uniform increase of this stress component 
throughout the sphere of the same size as a rr . The initial stress is thus isotropic. Note that crge shifts sign during the thermal 
equilibration process, in contrast to a rr . This can be understood in a physical picture: Consider the outer region that has been 
reached by the inflowing heat at a certain point in time. If that region were free it would expand, but it is kept in place by the 
inner unheated region that has not expanded thermally yet. This creates a negative stress on surfaces with normal at right angle 
to the radius vector. 

We conclude this section by a simple result. If one compares the instantaneous temperature drop Eq. (|26l and the instantaneous 
stress increase Eq. ( |3ll , one finds that the ratio is given by the adiabatic temperature — pressure coefficient: 



6T(r<R,t = 0) = (dT\ 



5p{r <R,t = 0) V d P J 



V/3s = -5- ■ ( 3? ) 



s 
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0.025 0.05 0.075 0.1 



t/x 

FIG. 3: The temperature of the sphere as a function of time for different radii normalized to the final temperature AT e9 = SQo / (Voc p ). After 
addition of heat at the surface, the temperature drops instantaneously throughout the sphere, showing adiabatic cooling. The time scale is given 
by the characteristic heat diffusion time, r — R 2 q/X . The longitudinal coupling constant is chosen to a; = 0.091. 




0.025 0.05 0.075 0.1 



t/x 

FIG. 4: The rr-component of the stress tensor, 1 — F2, as a function of time for a number of radii scaled with the initial stress, ao = 
I I§r Vbc"^ ' After the addition of heat at the surface, a rr immediately increases throughout the sphere. The stress is released as heat 
diffuses towards the center from the surface. 



B. The case when temperature is controlled at a mechanically free boundary 

The above studied case with heat-input control showed a rather simple cooling-by-heating behavior at short times or high 
frequencies. We now consider the case of controlling the temperature on the outer surface instead. There is still an effect, but 
it is not instantaneous. We only calculate the temperature in the center of the sphere. The surface is again mechanically free. 
Again, using the transfer matrix technique in the frequency domain, one finds that the temperature amplitude, ST(0, s) in the 
center is related to the temperature amplitude, 6T(R, s) at the surface by ST(0, s) = $>(s)6T(R, s), where 

x 3 — x 2 sinh(x) \ 
3a; [x cosh(x) — sinh(x)] — x 2 sinh(x) / 

Here x = st(oj), s = iu>. The characteristic diffusion time r(w) (Eq. d23l l) and thermomechanical coupling constant ai(ui) 
(Eq. (O) are in the general thermoviscoelastic case complex and frequency dependent and the temperature response can only be 
converted to the time domain numerically. 

In order to calculate the temperature signal as a function of time we again limit ourselves to the purely thermoelastic case, 
i.e., the case of a solid where r and a; are real and frequency independent. For a Heaviside temperature step at the surface of 
the sphere, ST(R, t) = ATH(t), the temperature at the sphere center is calculated via an inverse Laplace-Stieltjes transform of 
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0.025 0.05 0.075 0.1 

t/x 

FIG. 5: The #(9-component of the stress tensor, 1 — \F% + hFi, as a function of time scaled with the initial stress ao- There is an initial 
adiabatic positive step up in ogg throughout the sphere. The regions that have been reached by the incoming diffusive heat at a certain point in 
time experience a negative 00-stress component since the inner unheated regions pull the outer heated regions inwards. 



( GO 

ST(0, t)=AT\l-J2Rk cxp ( -v\- ) \ , (39) 



fe=0 




where the residues are given by 



Rk = 2(1 -cosj/fc) + y k siny k /(3ai) 

(1 - 3ai) cosy k + y k {2 - 3a ; ) siny fe ' 

Here the y k 's denote the positive roots of the transcendental equation 

cot(y) = - - JL . (41) 
V 3a/ 

In Fig. |6]we plot the solution Eq. d39l for various values for the coupling constant a; . Time is given in units of the characteristic 
diffusion time and AT = IK. We see that the cooling-by-heating effect is present also when a step in temperature (instead of 
heat) is applied to the surface. However, now the effect is not instantaneous but evolves gradually, reflecting the gradual heat 
diffusion at the surface mediated to the center by the stress field. Figure © furthermore shows that it is not enough to have 
a thermomechanical coupling (a p ^ 0) for the phenomenon to be present - only when q ^ c. p is there a cooling-by-heating 
effect. The next section studies the general, thermoviscoelastic case of frequency dependence of the response functions, which 
describes supercooled liquids. 
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FIG. 6: The temperature at the center of the sphere as a function of time. After the temperature at the surface has been raised, the temperature 
at the center of the sphere initially decreases. This only happens when q 7^ c p , i.e., when the longitudinal coupling a; is not negligible. The 
temperature step at the boundary is AT = IK, and the time scale is given by the characteristic diffusion time, r = R 2 ci/X. 



IV. THE THERMOVISCOELASTIC CASE 



The above time-domain results apply for a thermoelastic solid only, whereas the frequency-domain results are general. The 
thermoelastic examples handled so far only involved frequency-independent constitutive parameters, corresponding to the high- 
frequency (low-temperature) limiting values of the curves sketched in Fig. Q] However, Fig. Q] indicates that the value of the 
coupling constant a/ is larger at lower frequencies (higher temperature), thus suggesting that the effect of cooling by heating 
may be even larger in the very viscous liquid or simply at the glass transition. We investigate this issue in the time domain in 
this section. In the thermoelastic case the inversion of the problem to the time-domain could be made analytically. This is not 
possible in the thermoviscoelastic case where the constitutive parameters are complex frequency-dependent functions. 

In order to investigate the effect of going from solid to liquid we resort to numerical methods. Specifically, we transform 
$(s) of Eq. ( f38l > into the time-domain, accounting for the frequency-dependence of the constitutive parameters via t and a;. 
To do this we have to introduce a model of the constitutive parameters that enters via r and ai. It is common in rheology to 
illustrate models like the Maxwell model by rheological networks or even their electrical analogue. We use this approach to 
model the thermoviscoelastic behavior. The purpose of the model is to interpolate between the thermodynamic coefficients at 
high frequencies, kt,oo, c Pj00 , a Pt00 and at low frequencies kt.o, c p .o, ol p ,q- Network modeling assures internal consistency and 
agreement with the rules of linear irreversible thermodynamics. A one-parameter relaxation model implies the Prigogine-Defay 
ratio is unity, which is not the case for glucose. Rather, with T g = 300 K and Ac p = 1.14 • 10 6 JK -1 m~ 3 , Akt = 6.1 • 10 -11 
Pa - 1 and Aa p = 2.6 • 10~ 4 K _1 one finds 

„ Ac„Akt 

n = 5SL^ = 3 ' 4 ' <42) 

We are thus forced to consider a model with two relaxation elements that cannot be lumped into one. The model of Fig. [7] is 
suited for this purpose. In order to still make it simple, the two relaxing elements are taken to be Debye-like. The model has a 
simple mathematical formulation in the frequency domain. We change the independent variables compared to Eqs. (0 and © 
and consider the complex amplitude ST and Sp to be the controlled stimuli creating a linear response in the amplitudes 8 J and 
6e of the entropy density and dilation: 



Se) \ a p (uj) kt(u) ) \-Sp 



(43) 



Here ( p = c p /Tq, a p = Pv/Kt, Cp — Cv + Py/Kr and kt = 1/Kt- In the model the three measurable quantities, the 
isobaric specific heat, the isobaric expansivity, and the isothermal compressibility are related to the elements D, C, Ja(u), and 
JbM by 

C P = D 2 J a {lo) + C, (44) 
a p = -DJ A (u), (45) 
k t = J a (lu) + J b (lo). (46) 
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FIG. 7: Electrical equivalent diagram of the interaction of a volume element with its surroundings through two gates. The thermal gate where 
entropy displacement SS or temperature ST can be controlled, and the mechanical gate where volume displacement 5 V or pressure Sp can be 
controlled. SS and SV are generalized charge displacements and ST and — Sp are generalized voltages. The relaxational elements Ja and Jb 
are simple single relaxation time elements. 



The relaxing element J a is detennined by three parameters Ja,oo, AJa = Ja,o — Ja,oo and Ra- 



Ja = Ja, 



1 



i 

A. J A 



iujRa 



(47) 



and likewise for Jb- 

The parameters of the model can be established from the high and low frequency limits of £„, and a p . One finds D 



-AQ p /Aa p , J A ,o = -a Pt0 /D, J A ,, 
ratio of the model is given by 



ID, Jb.o = k t ,o - Ja,o and J S;OC = k T;OC - J A ,oo- The Prigogine-Defay 



n 



A( p An T 
{Aa p y 



AJ E 
AJa 



and the dynamic (frequency-dependent) Prigogine-Defay ratio [15] is given by 



A 



T" 
T" 



(48) 



(49) 



As expected, the Prigogine-Defay ratio is larger than unity, but becomes one if the element Jb is non relaxing, in which case 
the model reduces to a single-parameter model. In the case of J a being non-relaxing, the model degenerates with absence of 
relaxation in C p and a p . For glucose at 300K we calculate the parameters of the model to be D = — 1.46 • 10 7 PaK^ 1 , C ~ 
4.76-10 3 JK- 1 m- 3 , J A ,oo = 7.53- 10 -12 Pa -1 , AJa = 1.78 -1CT 11 Pa" 1 , J B ,oo = 8.55- lfT 11 Pa" 1 , and A J B = 4.33- KT 11 
Pa" 1 . 

The heat conductivity of glucose — needed to calculate the heat diffusion time — is A = 0.35 WK^'m -1 at 303 K according 
to Greene and Parks 111611 . The shear relaxation of glucose is modeled by a Maxwell model 



G(co) 



1 



i 

icurf 



(50) 



Here the high-frequency shear modulus can be taken to be Goo = 3.1 • 10 9 Pa JH]. The values of the thermodynamic parameters 
used to parametrize the model are given in Table U The temperature dependence of the shear viscosity causes the shift in the loss 
peaks of the relaxations. Parks et al. lfl7ll measured the viscosity of glucose in a wide temperature range from 295K to 418K. 
We fitted their tabulated data by the expression r](T) = 0.0125 exp((512.9K/T) 6 42 ) Pas. This holds within 20% over the entire 
temperature range except for the highest viscosity point of 9.1 • 10 12 Pas, which however is well beyond the glass transition and 
may be hard to measure reliably. A Vogel-Fulcher law is not as good a fit, deviating more than 40% in the measured temperature 
rang e. The rate parameters Ra = Ra(T) and Rb = Rb(T) are assumed to follow the temperature dependence of the viscosity 
111811 . It is found numerically that one should chose Ra(T) = Rb(T) = 30r](T) in order to get the loss-peak frequency of 
the shear modulus, G(u>), and isothermal bulk modulus, Kt{uj), to coincide. The relaxation of the different response functions 
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Quantity Value Reference 

oj ->■ 

k t ,o 15.4 ■ lCT^Pa -1 [4J 

c p , 3.05 • lf/jK^m- 3 [3J 

q p ,o 3.7-10- 4 K _1 [3] 

A (at 303K) O.SbWK" 1 ™- 1 \16} 

CO — i> oo 

k Ti oo 9.30 ■ lO^Pa^ 1 [4J 

Cp.oo 15.4 ■ lO^JK^m" 3 [3J 

a Pl00 1.1 ■ lO^K" 1 [3J 

Goo 3.1 • 10 9 Pa [5J 

TABLE I: Literature data for glucose (at 300K) used to parametrize the model depicted in the electrical equivalent diagram in Fig. [7] We have 
used a density of 1.52 • 10 3 kgm -3 to convert specific heat data from mass to volume. 



described by the model implies a frequency dependence of the longitudinal thermomechanical coupling via 

, f G(w)To<x 2 ,(u;) 

°'H = 7TT^ \ < V t v (51) 

(1 + ^G{w)k t {u:))c p {u) 
The modulus of this complex function was shown in Figure (fTJ. Also, the heat-diffusion time 

r(w) = R 2 (l-a l {u)))c p (w)/\ (52) 



now becomes complex. This makes the inversion to the time-domain non-trivial and thus Eq. (138b was inverted numerically. 
The algorithm for the inverse Laplace transform is an improved version of de Hoog's quotient difference method 11911 developed 
and implemented in Matlab by Hollenbeck f2(ill . 

The calculated temperature response in the middle of the sphere to a step of 1 K at the surface is shown in Fig. [8] Time is now 
scaled by the fixed real-valued diffusion time tq in the liquid regime, 

to = R 2 c pfi /X • (53) 

The figure shows that the effect of the thermomechanical coupling is absent at high temperatures. But as temperature is decreased 
and the liquid gets more and more viscous, a dip in temperature emerges. Going further down in temperature the phenomenon 
of cooling by heating becomes most pronounced slightly above T g . Even further down in temperature, in the glassy state, the 
effect is still present, but small. One may ask what happens if the expansivity a Pt00 vanishes, so that the phenomenon is absent 
in the glassy state: Will it still be present at the glass transition? The simulations shown in Fig. (O confirm this expectation. 
Putting a Pt00 = 0, but otherwise keeping the values of the rest of the parameters, we get a succession of temperature evolutions. 
Going down in temperature we see the cooling-by-heating phenomenon appearing at T g and afterward disappearing in the glassy 
state. Fig. [K)]shows the minimum temperature (STmin reached when a p>oc ^ as a function of temperature To, emphasizing the 
phenomenon as characteristic of the glass transition. 



V. EXPERIMENTAL VERIFICATION OF THE COOLING-BY-HEATING EFFECT 



To prove the existence of cooling by heating we molded glucose (a-D(+) glucose, 98%, Sigma- Aldrich) into spherical sam- 
ples with a thermistor placed in the center. Via the large negative temperature coefficient (NTC) thermistor we measured the 
temperature in the middle of the sphere. During measurements the glucose sphere is placed in a cryostat, that makes it possible 
to change the temperature at the surface of the sphere quickly compared to the characteristic heat diffusion time. A sketch of 
the setup together with a photo of one of the samples is shown in Fig. [TT] The photo shows the wires that lead into the sphere, 
connecting the thermistor to the terminals on the peek plate shown in the photo. When mounted on a holder, the terminals get 
connected to the multimeter that does the resistance measurement. 

The procedure in the experiments was the following. First we brought down the temperature to the desired starting level. Then 
we waited for the temperature to equilibrate. This was monitored by measuring the resistance every fifth minute. Typical waiting 
time was 18 hours. After the initial waiting time, we increased the sampling rate of the multimeter to about five data points 
per minute. This was done for one hour to get a baseline like the one in Fig. Q~2] Then we imposed a 5K temperature step and 
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0.04 



FIG. 8: The change in temperature, ST, at the center of a sphere, after a temperature step of AT = IK has been applied at the surface. A 
model of the thermoviscoelastic relaxation of glucose with realistic limiting thermodynamic parameters has been invoked. Time is scaled by 
the characteristic diffusion time To, which is 800 s for a glucose ball of radius 9.5 mm. The minimum of approximately — 5mK occurs for 
the 309K curve just above T g at 0.017i~o, corresponding to 14 s. Notice that although the cooling-by-heating phenomenon is present in the 
glassy solid state, it becomes more pronounced at the glass transition. The squares marking the minima at each temperature are plotted against 
temperature in Fig. QJj] 



K3 




0.04 



FIG. 9: A simulation similar to the one shown in Fig. ([§} except that a p ,oo has been put to 0, forcing cooling by heating to be absent in the 
glassy solid phase. It is seen that the cooling by heating effect still appears as a dynamic phenomenon at the glass transition. 



continued sampling data for another ten minutes with a sampling rate of 15 data points per minute. Fig. Q~2]shows the temperature 
measured by the NTC thermistor during a measurement with a step from 298K to 303K of the cryostat temperature. The baseline 
extends for about 20 minutes, then a characteristic temperature dip appears. The magnitude of the dip is 7.3 ± 0.2mK, which 
was reached 40 seconds after the temperature step was imposed. 

The experiment was repeated on three different samples; the inset in Fig. Q~2]shows the results of the first measurement done 
on each sample, at the same temperature. Each marker represents the lowest temperature reached in one measurement. They are 
plotted against time after the temperature step is initiated. It was not possible to reproduce the phenomenon on the same sample 
by recycling the temperature. We ascribe this to crystallization. Although the effect should be present also in the solid state, 
it is here considerably smaller and not observable with our temperature resolution. Nevertheless, the phenomenon was seen 
every time we repeated the experiment with a fresh supercooled sample. The sphere does not flow or deform to any appreciable 
degree even somewhat above the glass transition. The characteristic flow time ta ow is proportional to the viscosity and inversely 
proportional to the gravitational force rag. By a dimensional argument it follows that 

Tflow OC = T M , (54) 

pgr pgr 

which means that Tg ow ~ 10 7 Tm- Thus even at 310 K, where viscosity becomes 10 10 Pas and thereby the Maxwell time 3 s, the 
flow time is one year. 
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FIG. 10: The minima from Fig. [8]as a function of temperature. At high temperatures the phenomenon of cooling-by-heating is absent. Going 
down in temperature the liquid becomes more viscous and a dip in temperature emerges. Close to the glass transition the cooling-by-heating 
phenomenon gets most pronounced. As temperature is decreased further and the liquid enters the glassy state the effect is still present, however 
reduced in strength. 




0.4mm 



19mm 



FIG. 1 1 : Sketch of the experimental setup. The liquid is molded into a sphere, in which a small NTC-thermistor bead is placed in the center, 
connected to wires that lead to a multimeter that performs resistance measurements. The sphere is inserted into the cylindrical chamber of 
a cryostat. The photo shows one of the samples; at the time of the photo shoot (1 month after molding) the sample is no longer transparent 
because it crystallized. 



VI. DISCUSSION 

The concept of a longitudinal specific heat was identified in Ref. Qj as the relevant quantity within AC-calorimetric methods 
that utilize heat effusion. The principle of the simplest of these techniques ll2lll is to measure the complex temperature response 
T u at a plane surface to a heat current density jq^ generated at the same surface. The effusivity e = \f~\c is found from the 
measured specific thermal impedance Z = T^/jQ^ = 1/VwAc, and from the effusivity the specific heat can be calculated. A 
meticulous analysis of the thermomechanical equations of this problem showed that the specific heat that comes into play in this 
situation is q rather than c p . Effusivity measurements in spherical geometry have been shown also to involve the longitudinal 
specific heat iB E^Il . It is not a very well-known property, but it does appear in the textbook on elasticity by Landau and Lifshitz 
112311 . They show that the coupled thermoelastic equations decouple for certain boundary conditions of an infinite solid, namely 
when temperature at infinity is constant and deformation there is zero. They show further that the heat-diffusion equation is 
valid with a diffusion constant containing the effective specific heat ((1 + u)c p + 2((1 — 2er)cy)/(3(l — er)), where a is the 
isothermal Poisson ratio. Inserting a = (3Kt — 2G) / (6Kt + 2G) one readily finds that the effective specific heat is c\. The 
longitudinal specific heat also appeared in Biot's 1956 paper [7] in his diffusion equation for the entropy density. Although not 
very different from c p , there is a fundamental difference, and c; pops up in many thermoelastic problems when they are treated 
exactly. In particular, as we have seen in this paper, there is only a cooling-by-heating effect if c\ 7^ c p . We originally proposed 
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FIG. 12: The temperature in the center of a sphere of glucose, with diameter 19mm. At time t = 0s a step in temperature from 298K to 303K 
is made with the cryostat (see the setup in Fig, lilt. The temperature drops initially by 7.3mK as result of cooling-by-heating. The inset shows 
the minimum temperature STmim reached in three measurements, on three different samples, at the same temperature. 



the name longitudinal specific heat because this is the heat needed to increase temperature by 1 K if the associated expansion is 
confined to be longitudinal instead of isotropic. 

Transient thermal stresses induced by surface heating of a sphere have been considered theoretically by Cheung et al. ||24ll . 
Their interest was fragmentation of brittle solids by surface heating. A heat current was applied uniformly within the polar 
angle regime < 8 < 8q and the temperature and stress distributions calculated in time and space. This is a problem very 
similar to the one considered in this paper, although not generalized to situations with relaxation. However the cooling-by- 
heating phenomenon was not seen, since the standard decoupled heat-diffusion equation was used. The phenomenon thus seems 
apparently not to have been recognized in the literature, even though it belongs to classical continuum physics. 

Other kinds of phenomena have lately been termed by the phrase "cooling-by-heating" or related names. Thus in 1999 Aleshin 
et al. reported "heating through cooling" when a copper bar is first heated to 150 °C and subsequently rapidly cooled. This result 
in a temperature increase at the other end of the bar of about 4°C ll25ll . The authors presented the following explanation: The 
sudden cooling causes the bar to contract, producing an elastic wave that propagates towards the cold end of the bar. Under 
certain conditions this wave triggers a sequence of events responsible for an energy release, similar to a process called a "steam 
explosion" in an ionic liquid which is observed when a water jet interacts with a molten salt, an explosion that does not take 
place where the water hits the molten salt, but at the bottom of the container {26j| . In 2007 Zwickl et al. reported a "counter- 
intuitive cooling-by-heating" effect for laser cooling of a microcantilever, an observation they suggested is due to photothermal 
forces causing the lowest cantilever vibrational mode to cool while all other modes are heated |27ll. Very recently Man and 
Eisert discussed theoretically cooling of quantum systems by means of incoherent thermal light Il28ll . While coherent driving of 
a quantum system can mimic the effect of a cold thermal bath, the novel idea is that under certain conditions even incoherent 
"thermal" light can be used to cool a quantum system. 



VII. SUMMARY 



We have shown that cooling by heating occurs at the center of a solid spherical sample if it is heated at a mechanically 
free surface, reflecting a non-trivial thermomechanical coupling where the temperature initially decreases in the interior of the 
sphere. What happens is that, as heat diffuses into the outermost parts of the sphere, these parts expand and build up a negative 
pressure at the center of the sphere. This negative pressure couples to the temperature via the adiabatic pressure coefficient 

(l5p ) ■ T ne PP 0S it e effect also applies, of course: if the temperature of the surface is lowered, heating by cooling will be 

observed. In ordinary solids the cooling-by-heating effect is almost negligible because their thermal expansion is generally 
small. The effect is particularly large for liquids close to their glass transition. The cooling-by-heating phenomenon establishes 
the difference between the longitudinal and isobaric specific heat, since the effect is only present when these two quantities differ. 
This is the case when the shear modulus is non-vanishing compared to the bulk modulus and, simultaneously, the isobaric and 
isochoric specific heats differ significantly. Analytical results show that the phenomenon occurs in the elastic case (the solid), 
and numerical results based on a model of the glass transition with parameters determined by glucose data show that the effect 
is present dynamically also in the very viscous liquid. Even in the hypothetical case when the glassy state is assumed to have 
zero expansivity, the phenomenon still appears at the glass transition. 
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The numerical results based on glucose data indicate that the drop in temperature at the center of the sphere is of order 5 mK 
with a duration of approximately 15 seconds when the temperature is increased by IK on the surface of a sphere of radius r = 10 
mm. This prediction was confirmed experimentally. 
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Appendix 

1. Solution from section Hill - the transfer matrix method 

The solution to the problem given in section [HI] is based on the transfer matrix formulation [0 |H |29| of the general solution 
of the thermoviscoelastic problem in a spherically symmetric case. Including the radial stress field and the time-integrated 
heat-current density that relate to the temperature and displacement fields, the authors of Refs. [J, end up with four coupled 
equations to solve. Laplace transforming the equations relating the four fields and solving the resulting inhomogeneous system 
of four ordinary differential equations, the result is in the general form of a transfer matrix T(j, i) that links the dimensionless 
complex amplitudes of the fields at the boundary r, with those at : 



/ 8Pr\ 

ST 
SV 
\SSj 



/ 5p r \ 
Sf 
SV 

\ssj 



(A55) 



Here SS, SV, ST, and Sp r are the complex amplitudes of entropy, volume, temperature, and the radial component of pressure 
(5p r = — oy r ), respectively. The elements of the transfer matrix are given in reference $2^. From this general solution one can 
work out different cases, like the ones in Sec. [HI] The boundary condition at f i = 0, giving net flux of heat through the center of 
the sphere, was set equal to zero, SSi = 0, since heat is supplied uniformly across the surface f 3 giving a spherically symmetric 
case. For the same reason, SV\ = 0. At the mechanically free outer boundary f 3 the heat supplied, SS 3 , is given and 5p r3 = 0. 
Letting f = f 2 be an intermediate variable radius between f 1 and ¥■$ one has 



(5p r \ 

sf 

SV 
\S~SJ 



= f(2,l) 



fSp r \ 

sf 





(A56) 



Also 



ST 
SV 
\SSJ 



T(3,l) 



/Sp r \ 
Sf 


V y , 



(A57) 



This leads to Sp r> \ = 

^2(3,1)^(3,1) u ~ 

This implies 



Ti2(3,l) 
T11 (3,1) 



<m, or Sf 1 - -pSmSprt whereby SS 3 = (f 42 (3, 1) 



T , 4i(3,l)Ti 2 (3,l) 
Tii (3,1) 



)tffi = (T 4 i(3,l)- 



(3,1) 



T 41 (3,1)T 12 (3,1)-T 42 (3,1)T U (3,1 



■)SS 3 



(A58) 



Sp r 



Ti 2 (3,l) 



T 4 i(3, 1)T 12 (3, 1) - T 42 (3, l)r n (3, 1) 



)SS 3 



(A59) 
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Inserting this into Eq. (IA561 > one has 

5p r (f) = Sp r ^2 



f 11 (2,l)f 12 (3,l)-f 12 (2,l)f 11 (3,l) 
r 4 i(3, l)f 12 (3, 1) - T 42 (3, l)T n (3, 1) 



)6S 3 



(A60) 



Sf(r) = Sf 2 = ^(M)T 12 (3,l)-f 22 (2,l)f ll( M) 
~ T 41 (3, 1)T 12 (3, 1) - T 42 (3, 1)T U (3, 1) 



(A61) 



(r) = ^ = ^(2,l)T 12 ( 3 ,l)-T3 2 ( 2 ,l)T ll(3 ,l) 
T 4 i(3, 1)T 12 (3, 1) - T 42 (3, 1)T U (3, 1) 



(A62) 



<J5(f ) = 55 2 



Inserting the actual explicit values of the transfer matrix elements from 
f = r 2 R = f3, yields 



f 4 i(2, l)f 12 (3, 1) - f 42 (2, l)f n (3, 1) - 
T 4 i(3, 1)T 12 (3, 1) - T 42 (3, l)fu(3, 1) " 

and evaluating them in the limit of f\ 



Sp r (f) 
ST(r) 

SV(f) 
(55(f) 



35g 
5(1+5) 
1 

£ 



fcosh(f) — sinh(f) 



1 

R 3 r 3 (R cosh(R) - sinh(ij)) 
3a 2 g sinh(f) 



SS(R) 



R3 & 2 g + d{l+g) 



3 ~,~2 



SS(R) 



c{l + g) 



f{R cosh(R) - sinh(i?))_ 
c(l + g)) f cosh(f) — sinh(r 



& 2 g + c(l+g) R cosh( J R) - sinh( J R) 



5S(R) 



fcosh(f) — sinh(f) 



Rcosh(R) - sinh(i?) 



^SS(R) 



(A63) 
0, putting 

(A64) 
(A65) 
(A66) 
(A67) 



The transformation back to dimensionalized physical quantities is performed by noticing that f = kr,R = kR, u = ku, 
where k = y/sqj\. Furthermore c = T ci/K T , a = T a p , g = 4G/(3A" T ), a rr = <r rr /K T , SV = 5Vk 3 /{4ir), SS = 
5Sk 3 1 (AttKt), ST = ST /To. Since the entropy-displacement is positive in the direction of r, it is related to the heat input at 



the outer surface by SQ = —TqSS, opposite of the convention in Ref. 0- Recalling that 
(1241 and <|25) from Eqs. (|A65T > and (lAlvgt , 



f V one now easily derives Eqs. 



2. Inverse Laplace-Stieltjes transforms. 

If a stimulus <fi s e st {s = iuS) on a linear system gives rise to a response 7 s e st , where 7 S = f(s)<fi s , the response to a Heaviside 
input 4>oH{t) will be -y(t) = F(t)(/)Q, where F(t) is the inverse Laplace transform of f(s)/s (or the inverse Laplace-Stieltjes 
transform of /(s)). We perform the inverse Laplace transform via the calculus of residues as 

F(t)= V Res(^,s n )e 5 " t (A68) 
^— ' s 

poles s n 

Put p = r/R and t = (c;/A)i? 2 . Then kr = y/srp. Choose furthermore initially time units such that r = 1. Then the two 
expressions of Eqs. (l2TT i and (l22l when divided by s become 

fi(s) _ smh(p^s) VjS 

s Sp\fs y/s cosh(-y/s) — sinh(-^/s) 

/2(a) 1 p^/s cosh(pyi) - sinh(py / 5) 

s p 3 s \fs cosh( v / s) — sinh(y / i) 

Both expressions have simple poles at s = with residue 1. The other poles are on the negative real 3.X1S S n — X, r , where 
x n are the positive roots of the transcendental equation tan(x) = x. x\ ~ 4.493409457 and a; 2 sa 7.725251836 and x n w 



(A69) 
(A70) 
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\J {it/2 + nn) 2 — 2 better than 1 ppm for n > 3 . The residues now become, respectively 




(A72) 



(A71) 



and the corresponding time-domain functions 



Fi{p,t) 




sa\(px n ) 
sin(x„) 




(A73) 



F 2 (p,t) 



1 + 2 




sm(px„) - px n cos(px n ) 



x^ sin(x„) 



(A74) 



Using these expressions and reintroducing the characteristic heat diffusion time r, one derives Eqs. (l32l . d33l l. and (l34l l. 
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